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Electrostatic polarization is important in many nano-/micro-scale physical systems such as col- 
loidal suspensions, biopolymers, and nanomaterials assembly. The calculation of polarization po- 
tential requires an efficient algorithm for solving 3D Poisson's equation. We have developed a useful 
image charge method to rapid evaluation of the Green's function of the Poisson's equation in the 
presence of spherical dielectric discontinuities. This paper presents an extensive study of this method 
by giving an convergence analysis and developing a coarse-graining algorithm. The use of the coarse 
graining could reduce the number of image charges to around a dozen, by 1-2 orders of magnitude. 
We use the algorithm to investigate the interaction force between likely charged spheres in different 
dielectric environments. We find the size and charge asymmetry leads to an attraction between 
like charges, in agreement with existing results. Furthermore, we study three-body interaction and 
find in the presence of an external interface, the interaction force depends on the curvature of the 
interface, and behavior a non-monotonic electrostatic force. 



I. INTRODUCTION 

Electrostatic interaction is of fundamental importance 
in many fields of sciences at nano/ micro scales. The long- 
range and many-body characters of the interaction lead 
to various challenges in the understanding of complex 
electrostatic phenomena One debate in electro- 

statics is the like-charge attraction (LCA) in colloidal 
science, which is not in agreement with the prediction 
based on the classical Poisson-Boltzmann (PB) theory. 
The LCA has been observed in a lot of experimental stud- 
ies motivating wide computational and theoretical 
interest to understand Coulomb many-body phenomena 
(see @, 0, HI for reviews). 

One of the central issues is the polarization effect due 
to the presence of dielectric interfaces. The dielectric con- 
stant of nano/microparticle materials could range from 
about 2 (hydrocabon objects) to the infinity (metallic 
objects), while that of surrounding environments could 
range from 1^2 (air) to 80 (water). The pairwise in- 
teraction potential and the ion distribution surrounding 
nanoparticlcs significantly depend on the dielectric prop- 
erties of the system. The effect of dielectric disconti- 
nuities has been often studied in various systems such 
as colloidal suspensions, liquid droplets in clouds, self- 
assembly of nanoparticles, and biopolymers. 

In computer simulations, the polarization effect due to 
dielectric interfaces remains a technical difficulty, and the 
algorithm development remains up-to-date attention 0- 
Il2| . Direct numerical solutions of the Poisson's equation 
with finite element methods [HI are intensive for Monte 
Carlo and molecular dynamics simulations( for a survey 
of electrostatic algorithm in simulations, see [3), and 



analytical solutions are only available for simple geome- 
tries such as one spherical interface for which the method 
of image charges can be applied (l5l - [i"8| . 

The electrostatic interaction between multiple dielec- 
tric spheres usually needs to use bi-spherical harmonic 
expansion [i~9l - [2l1 | . Spherical harmonics are computation- 
ally expensive and slowly convergent for ions approaching 
to the surfaces in simulations |22l). Other methods have 
also been developed [l2l.l23l.l24l]. But, mostly analytical 
studies limit to no more than two spheres. Alternatively, 
we have extended the discrete image approximation [25[ 
of the Neumann's image principle to multiple spheres by 
reflections of discrete image charges among different in- 
terfaces [26j. In this paper, we will present an extensive 
study and aspects such as the coarse-graining strategy 
and the convergence analysis will be proposed. Further- 
more, we investigate the interaction force between likely 
charged spheres under size and charge asymmetries, and 
in the presence of an external interface, where the sur- 
face charge distribution is modeled by discrete charges. 
The phenomenon of LCA at the conducting limit was 
widely known, and this work provides an efficient tool 
for a quantitative study of many-sphere interactions at 
finite dielectrics. 



II. METHOD 

We are interested in a system made of a cluster of 
dielectric spheres and a large number of ions with space 
charge distribution p(r) = J^j Qj^( r ~ r j)i where 5 is the 
Dirac delta, r 3 and qj is the location and charge of the 
jth ion. The electric potential $ of the system satisfies 
the Poisson's equation, 



'Electronic address: xuzl@sjtu.edu.cn 



- V • e(r)V$(r) = 47rp(r). 
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We suppose the permittivity inside the dielectric spheres 
is £j which could be different for different spheres, and in 
the bulk solvent it is e , then e(r) is piccewise constant. 
Eq. ([1]) implies the boundary conditions across the in- 
terfaces, namely, the electric potential $ and the electric 
displacement £(r)<9$(r)/9n are continuous on the inter- 
faces. 

Suppose the electric potential decays to zero at the 
bulk solvent far way from the dielectric spheres. The 
solution of Eq. ([1} can be written as [27} , 

$(r) = [ G{v,v')p{v')dv' =5^(1-,^), (2) 

J 3 

where G(-, •) is called the Green's function of the Pois- 
son's equation, described by, 



V • e(r)VG(r, r') = 4vr(5(r - r'). 



(3) 



Physically, the Green's function is the electric potential 
produced by a unit point charge at the source point r'. 
Once the Green's function is explicitly given, we ob- 
tain the solution of the Poisson's equation for an arbi- 
trary charge distribution. As the source charges are lo- 
cated outside spheres or on the interfaces, we rewrite the 
Green's function outside the dielectric spheres as the sum 
of the Coulomb potential and the polarization potential, 



G(r, r') = $ cou i(r, r') + $ pol (r, r'), 



(4) 



with $ CO ui(r) = l/£o|r — r'[, the Green's function in free 
space. 

A. Method of image charges 



from the origin to the point of the Kelvin image along 
the radial direction. Using a Gauss-Legendre quadrature 
to the continuous line integral leads us to discrete image 
charges for approximating the polarization potential [25( , 
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where q\ = —^R/r 1 and Xi = vk are the parame- 
ters of the Kelvin image charge, and q m = uj m ^R/(2r') 
and x m = rx {(1 — s m )/2} 1 ^ cr , for m = 2, ■ • ■ ,M, are 
the charges and positions of the discrete point images 
due to the line integral contribution. The coefficients 
{w m , s m , m = 2, 3, • ■ ■ , M} are the (M — l)-point Gauss 
weights and locations on the interval [—1,1], which can be 
found in many literature; for example, Numerical Recipes 
[28| . Besides the Gauss quadrature, the integral can be 
approximated by a hypergeometric function [111 ] , and his- 
torically other approximate methods of images p9l - [3l| 
are also developed. 




For the problem of a unique dielectric sphere in the 
system, it is know the Green's function solution can be 
obtained in the spherical coordinates or by the image 
charges, and our discussion focus on the later one. Sup- 
pose R is the radius of the sphere, r and r' are radial 
distances of r and r', and the spherical center is the co- 
ordinate origin. The Neumann's image principle (l6j ex- 
presses the polarization potential as a line integral (See 
also [l3, for a recent discussion) , 

1 f K f(r K /x) j 

with 

f(t) = - 1 6(t-l) + ^H(l-t). (6) 

Heree = £i/e , 7 = (£i-£ )/(£i+£o) and a = £ /(£i+£ ) 
are three real constants, and rx = rxr'/r', tk = R 2 /r', 
and x = xr'/r'. It can be seen the Dirac delta function 
in the charge density corresponds to a point charge called 
the Kelvin image at rx, and the Heaviside function H(-) 
corresponds to a continuous line charge. The integral is 



FIG. 1: Illustration of the construction of image charges by 
reflections among surfaces. A source charge q produces the 
lst-level image charges inside each dielectric sphere, symbol- 
led by blue circles and one black circle. Each image produces 
the 2nd-level image charges inside other spheres, symbolled by 
red rectangles for the black lst-level image (Only this image's 
2nd- level image charges are shown). The series is convergent 
by recursively performing this procedure. 

When there are multiple dielectric interfaces in the sys- 
tem, the polarization potential is constructed by reflec- 
tions between different surfaces, as is illustrated in Figure 
Q] Suppose we have Nd dielectric spheres. First of all, for 
a source charge, M image point charges are generated in 
each dielectric sphere by using the image expression for 
one sphere, Eq. ([?])• The potential due to the summation 
of the source charge and the image charges in a specified 
sphere satisfies the interface conditions on this spherical 
surface, but fails to satisfy the interface conditions on 
the other spherical surfaces, and therefore M new im- 
age charges in each of the other sphere have to be used 
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for each image charge. This procedure is recursively per- for the polarization potential of the whole system: 
formed (see Figure[l}, leading to the following expression 
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where the subscript "jkmkjk-imk-i • • • jimi" refers to 
the next-level image charges of the charge with subscript 

"jfc-l<7Jfc-l • • • JlTOl". 

In Eq. (JS]) , there are NdM image charges in the first 
term, which is double sums and we call they are the im- 
ages at the first level. The nth term is called the nthe 
level, which is 2n-fold sums and with N d (N d - l)™" 1 ]^™ 
images. We truncate the series at Lth level, $ po i « 
$1 + $2 + • • • + where <6j is the 2v>fold sums for 
images at the ith level, then the total number of image 
charges is, 



N d M L+1 {N d - l) L - N d M 
N d M - M - 1 ' 



(9) 



A simple truncation of the series may lead to low ac- 
curacy and expensive computation. It could be noticed 
that, in the one-sphere problem the Kevin image charge 
has an opposite sign of the image line density, and overall 
they are electric-neutral and behavior as a dipolc. Higher 
levels give higher order poles, and presumedly the series 
((HJ has a very fast convergence. 

All image charges are concentrated inside spheres with 
a high density. It favors to combine some of them to- 
gether by a coarse-graining strategy. This will greatly 
reduce the number of image charges, as is discussed later 
on. 



B. Convergence 

We consider a special case that two spheres of the same 
radius R with separation d = ri2 — 2i?, and ri2 is the 
distance between two centers. The convergence rate as 
function of e — d/R is estimated. We calculate the self 
energy of a unit source charge, which is equivalent to 
calculating the polarization potential at the source point. 
The source charge is placed at the middle of two spheres, 
and the convergence is faster if the source charge deviates 
from this point. For asymmetric spheres, the following 
estimation still holds by taking e = d/ R> where R > is the 
larger radius. Therefore we need only do the analysis for 
the special case to get a lower bound of the convergence 
rate. 

We take M = 2 so each reflection of a charge results 
in an image dipolc. Suppose an image charge qi at n 



distance to the center of the sphere including itself. Let 
r s = ri2— r\, then the potential contribution of this image 
charge at the source site is, 
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(10) 



Then this image's image dipole within the other sphere 
have the potential contribution, 



-jRqi/e Q 
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(11) 

Let us define y = r s /R, then recalling e = d/R we have, 
e 2 + Ae + 2 



2 + e 



< V < 2 + e. 



(12) 



We calculate the absolute of the ratio between two po- 
tentials = $i/$n, which is 

M (i-f) -(y- 1 -!) 

By using inequalities (|12|) . a simple derivation gives us 
an estimation of the convergence rate, 



m > vV2(v-mv-v) + v(2y-i~§) ( . 

1 u W(i-f) [ ] 

where r) — 1 + e/2. 

At the limit of infinite radius, e = and ry = 1 and the 
convergence rate is bigger than a constant value, 1/|7|- 
And at a certain separation distance, for example, e > 1, 
the convergence of the image approximation is fast. 



C. The coarse-graining algorithm 

After a couple of reflections, mostly image charges will 
be within the half radius of each sphere. We develop a 
coarse-graining technique to reduce the number of image 
charges, which is useful for dynamical simulations. The 
algorithm is as the following. 

Let (5 > be a given error criteria. Suppose i"j. m and 
Tj <n are the location of two image charges due to the 
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source j within sphere k 7 which has center Ofc and radius 
. If they both have radial distances less than half the 
radius, i.e., \rj iTn - O k \ and |r,-,„ - O k \ < 0.5i? fc , we 
calculate the distance between two points 



(15) 



This distance is compared with the product of the dis- 
tance criteria and the radius, and if d mn < RkS, then we 
combine these two charges together, and place the new 
image charge at the weighted center of the two charges. 
Let w m = \qj >m \/(\qj,m\ + \lj,n\) and w n = 1 - w m , then 
the charge position is at w m rj, m + w n Tj in , and the charge 
strength is trivially the sum of two charges, qj m + qj. n - 
The performance of the algorithm will be discussed in 
the following section. 



III. RESULTS 

In this section, the interaction between charged spheres 
with different dielectric permittivities is investigated. We 
suppose the surrounding environment of the spheres is 
a homogeneous dielectric, which can be considered as a 
mimic of highly dilute solutions or a dielectric material 
in vacuum. 

We consider two spheres of radii R\ and R2 with scaled 
separation d — (ri2 — Ri — -R2V-R1, each has N = 60 
charges on its surface. The electrostatic energy of the 
system is composed of source-source interactions and the 
source-image interactions, 



- 1 V r ^ k 



<i,<i>j- 

•i-rj, k \ 



(16) 



where k ranges over all images of the jth source charge, 
and the image strengths and locations depend on the 
separation distance. The prime at the northeast of the 
summation represents that when i = j the source-source 
interaction and the interaction between source and its 
first few images which do not change with d should be 
not included. We define the interaction force between 
spheres by the negative gradient of the energy along the 
separation, 



F(d) = -dU(d)/dd. 



(17) 



Positive value of F represents a repulsion between two 
particles, and negative one represents an attraction. We 
use dimcnsionless unit by setting e Q = 1, the radius i?i = 
1 and interfacial ion valence qi = Qi = 1/N for the 
first sphere. The system asymmetries are described by 
e = £i/s: , R = R2/R1 and Q = Q 2 /Qi- 

The configuration of interfacial ions are generated from 
a Monte Carlo simulation of Coulomb repulsive particles 
on the surface. With the zero-temperature limit, the par- 
ticles form a quasi Wigner-crystal structure, which is con- 
sidered the ions are uniformly embedded on the surface 
in our study. 



A. Validation of accuracy for different parameters 

We let R = 1, Q = 1 and two sets of dielectric ratios 
e = 0.02 and 20. These sets of systems are used to test 
the convergence and accuracy of the approximation with 
image charges. We compare the solutions between dif- 
ferent numbers of image charges for one-sphere problem, 
M, and deeps of reflection levels, L. 
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FIG. 2: Convergence of the interaction force between two 
charged spheres as a function of separation. Solutions with 
large L are used for the "exact" solutions. R = 1, and Q = 1. 
(a) e = 0.02; (b) e = 20. 

The results are illustrated in Figure [2J where the "ex- 
act" is the reference solution from large L and M. It is 
seen all curves are very accurate for a separation d > 0.5. 
For small d, however, more images should be used. The 
approximation has a fast convergence to M, and M = 2 
is already a useful accuracy since the charge neutrality 
is satisfied. In regard to the parameter L, the conver- 
gence has a dependence on e. For small e, the image 
solution converges rapidly with L. For large e, an odd 
L gives a fast convergence. This makes sense because 
the Kelvin image shares the same sign as the source for 
e < 1. While for e > 1, the Kelvin image has an opposite 
sign, requiring the truncation at odd L. 

To test the performance of the coarse-graining algo- 
rithm, we calculate the same system as in Figure [2] for 
the case e = 20, and two sets of parameters L = 3, M = 3 
and L = 5, M = 3, denoted by L3M3 and L5M3, which 
lead to 78 and 726 image numbers, respectively. We take 
S = 0.001, 0.005 and 0.01. The number of image charges 
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FIG. 3: Mean number of image charges of each source with 
the coarse-graining algorithm. The image numbers for L3M3 
and L5M3 without the coarse graining are 78 and 726, re- 
spectively. The inset plot shows the relative errors of the 
L5M3 forces with the coarse graining in comparison to those 
of original image charges. 



and the relative errors of the force deviated from the 
original image charges as a function of the separation d 
is illustrated in Figure [3] The coarse graining reduces 
the numbers of L3M3 and L5M3 to around 15 and 90, 
decreased by a factor of 5 and 8, respectively, in the case 
of d < 0.4. The image number is further decreased for 
a larger d, as is shown that most of curves drop down 
to less than 10 with the increase of d. For the relative 
errors of the forces of L5M3 with the coarse graining 
compared to those of original image charges, all three 
S gives high accuracy, generally, less than 1%, except a 
point at d = 0.725 for 8 = 0.01. 

In the following calculations, L5A/3 will be used, for 
which the image method has a reasonable accuracy for 
both small and large dielectric contrast. 



B. Size and charge asymmetries 

Many literature has documented that likely charged 
spheres attract one each other due to the polarization 
under different size and charge, and different dielectric 
contrast. In a recent work, Lekner (32l | demonstrated 
that at small separation two conducting spheres almost 
always attract each other. For finite dielectrics, Bichout- 
skaia et al. [l2[ calculated the dependence of interaction 
force on size, charge, and dielectric ratios; see also ref- 
erences therein for an overview of historic works. Most 
of existing calculations use uniform representation of the 



surface charge, instead of the discrete representation in 
our work. 

We calculate the force for varying dielectric ratios by 
taking two groups of settings: one is with R = 1 and Q = 
3, and the other one with R = 3 and Q = 1. The results 
with size and charge asymmetries for dielectric contrast 
e varying from 0.02 to 1000 are shown in Figure |U The 
plots demonstrate the breaking of symmetry does lead 
to an attraction between likely charged spheres at small 
separation. When e > 25, the dielectric spheres have very 
similar property as conductor. In the right panel of the 
results, the attraction regime for high dielectric contrast 
is about R, which is because the second sphere has a 
radius of 3R and thus the image dipole in this sphere has 
longer-ranged interaction. 



C. Three-body interaction 

The polarization effect is essential a many-body in- 
teraction. The presence of the third particle affects the 
interaction force between two charged spheres. We study 
the interaction of two charged spheres with R = 1 and 
Q = 1 in the presence of an interface of different curva- 
tures, i.e., another dielectric sphere with varying radius. 
The external sphere has a factor of Re of the size R\, 
so there are three dimensionless distances, d, R and Re- 
The sphere is placed at at a distance R + Re away from 
the middle of the charged spheres; See Figure [5] (a). The 
external sphere is charge-neutral, but with the same di- 
electric contrast to the background, e, as the charged 
spheres. As the interactions with low dielectric ratios 
are always predicted repulsive, we focus on a dielectric 
permittivity e = 80, which models water droplets in the 
air H|. 

From Figure[5](b), it is observed the properties of these 
two spheres at a specified separation strongly and non- 
monotonously depend on the size of the external inter- 
face. Interestingly, all curves show attraction at small 
separation. When the separation is at the range of spher- 
ical radius, both large and small Re lead to repulsive 
forces (Re — 2 and Re > 50), and the intermediate radii 
(Re = 5 and 10) are attractive. Starting from d = 0.5, 
most of curves predict attractive forces for a wide range 
of d; but for small Re = 2 they are always repulsive. It 
could be interesting that, at Re = 50 and 500, the two 
particles have firstly an attraction and then an obvious 
repulsive force, and then they are attractive which con- 
tinues for a long range. Finally, all cases are repulsive at 
the far field (not shown in the figure) as the polarization 
force is essentially a dipole interaction, shorter than the 
direct Coulomb repulsion. At the planar limit of the in- 
terface, Re = 10 5 , the charged spheres are always repul- 
sive (except the short separation since the surface charge 
distribution is not completely symmetric), in agreement 
with existing study [Hj]. 

We use this model to validate the coarse-graining al- 
gorithm again. In Figure [5] (c) and (d), the results are 
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the image number and interaction forces for the method 
with the coarse graining with the parameter S = 0.005, 
and we see Figure [5] (d) accurately reproduce the re- 
sults without using the coarse graining. For problems 
with three dielectric spheres, the total number of image 
charges with L5A/3 for one source charge is 13995 by Eq. 
(|9"|) . Shown in Figure [5] (c) , the numbers of images are 
generally around 1550 for Re > 5 and when the separa- 
tion is not extremely small. The reduction is significant, 
but it is still high because most of image charges in the 
large sphere can not be coarse-grained. If the size of the 
external sphere is not so large, for example, in the case of 
Re = 2, the number becomes less than 30 when d > 4.5. 

IV. CONCLUSIONS 

In summary, we have developed a fast algorithm for 
evaluating electrostatic interaction between multiple di- 



electric spheres and reported the results for pairwise in- 
teractions of charged spheres in different dielectric envi- 
ronments. It is shown the algorithm is attractive for ac- 
curate approximation of the polarization potential. And 
it is found that the like charge attraction is a general phe- 
nomenon for a wide range of dielectrics in condition of 
breaking symmetries such as size and charge asymmetry, 
and in the presence of an external interface. 
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FIG. 5: Electrostatic interaction between like charged spheres near a neutral interface, e = 80 for all three objects (two charged 
spheres and a neutral large sphere), (a) Schematic setting; (b) Interaction force between two charged spheres along the radius as 
a function of separation d with L5M3 for image charges, (c) Mean number of image charges with the coarse-graining algorithm 
of parameter 8 — 0.005. (d) The interaction force with the coarse-graining treatment for image charges. 
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